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ABSTRACT 

Ereely decaying relativistic force-free turbulence is studied for the first time. We initiate the magnetic 
field at a short wavelength and simulate its relaxation toward equilibrium on two and three dimensional 
periodic domains, in both helical and non-helical settings. Eorce-free turbulent relaxation is found to 
exhibit an inverse cascade in all settings, and in 3D to have a magnetic energy spectrum consistent with 
the Kolmogorov 5/3 power law. 3D relaxations also obey the Taylor hypothesis; they settle promptly 
into the lowest energy configuration allowed by conservation of the total magnetic helicity. But in 2D, 
the relaxed state is a force-free equilibrium whose energy greatly exceeds the Taylor minimum, and 
which contains persistent force-free current layers and isolated flux tubes. We explain this behavior 
in terms of additional topological invariants that exist only in two dimensions, namely the helicity 
enclosed within each level surface of the magnetic potential function. The speed and completeness of 
turbulent magnetic free energy discharge could help account for rapidly variable gamma-ray emission 
from the Crab Nebula, gamma-ray bursts, blazars, and radio galaxies. 

Subject headings: magnetohydrodynamics — turbulence — magnetic fields — gamma-rays: bursts — 


1. INTRODUCTION 

The most extreme sources of high energy astrophysi- 
cal radiation are widely believed to exist in magnetically 
dominated, relativistic environments. Jets powered by 
super-massive black holes, plasma winds driven by pul¬ 
sars, and gamma-ray bursts are prime examples. The 
violent intermittency of gamma-ray production by these 
systems could be taken as strong evidence that turbu¬ 
lence is critically linked to their radiative output. And 
yet, the physics of magnetically dominated, relativistic 
turbulence remains nearly unexplored. 

The importance of understanding turbulence in this 
new regime is underscored by the discovery of power¬ 
ful gamma-ray fiyes originating within th e Crab Neb- 


Abdo et al. 2011 Tavani et al. 2011). Moreover, 


time variability seems to be ubiquitous a mong 


gamma-ray emitte rs; the bla zars PKS 2155-3 04 (Aha- 
roni an et al.[|2007|), 1510-0 89 (jSaito et al.||2013|), and 3C 
279 (Ha yashida et al.|2015), as well as radio galaxies such 
as M87 ( [Aharonian et al.|20 06) and IC 310 (Aleksic et al. 


2014), have each been observed to produce sporadic, high 


intensity outbursts of gamma-radiation. Such dramatic 
enhancements of synchrotron or inverse Compton emis- 
sivity require a reservoir of free energy to spontaneously 
energize the active region’s electron population. If that 
free energy resides in magnetic fields, then its discharge 
could be triggered by magnetic reconnection — the gen¬ 
eral pi cture of which has been rendered in many diflFerent 
ways (iLyntikov fc UzdensT^ 2QQ3|^ Lazarian et^ al.| 2003 


Zhang ^ Yan||2Ulll IMcKinney Uzdensky||2U12[ |Sironi| 

& Spitkovsky 120141 |Blandford et al.||2015|). 

In this paper we intend to demonstrate that magnetic 
free energy discharge can proceed from field geometries 
that are far more general than those typically consid¬ 
ered in reconnection models, and on a time scale that 
is not limited b y the rate with which micro physical or 
anomalous (e.g. Lazarian & Vishniac 1999) resistivity 
can destroy magnetic flux. T'his amounts to extend- 


ing the historical problem of magnetic relaxation (e.g. 
Chandrasekhar & Eermi 1953) to relativistic, magneti- 
cally dominated conditions! We focus on only a few of 
the many aspects of this topic that could be studied. 
Briefly, they are: (1) the rate and completeness of mag¬ 
netic free energy discharge in various topological settings, 
(2) a characterization of persistent non-linear structures, 
and (3) the spectral energy distribution of freely decaying 
relativistic force-free turbulence. To be most relevant for 
astrophysical gamma-ray emission, we are interested in 
regions far from any solid boundaries that could anchor 
the magnetic field (so periodic domains are appropriate), 
and where the plasma is nearly perfectly conducting, in- 
viscid, and magnetically dominated — conditions which 
are the domain of force-free electrodynamics (EEE) the¬ 
ory. 

Eorce-free electrodynamics forms the basis for histori- 


lianj 

IT^ 

Spitkovskyj 2006 

) and angular momentum ex 

traction from black holes 

(Blandford & Znajek 1977) 


ing these highly relativistic settings 


2010 Yang et al. 2015 Gralla et al 


Palenzuela et al. 
^OT^. It can be 


derived from relativistic magnetohydrodynamics (MHD) 
when the electromagnetic contribution to the stress- 
energy tensor greatly exceeds contributions from matter, 
and hence it captures the essential non-linear dynamics 
of relativistic MHD for the regime of interest. It also 
admits a numerical approach that is more robust and 
efficient than relativistic MHD solution schemes. 

Turbulence in force-free electrodynamics has only been 
considered in a few previous studies. The theory of 
Alfven wave turbulence in the presence of a strong guide 


field, originally formula ted for Newtonian MHD by |Gol-| 
dreich & Sridhar| (1995), has been extende d to the mag¬ 
netica l ly dom inated, relativistic regime by [Thompson fc:| 


Blaes ([1998 ). Alfven wave turbulence has since been 


studied nunierically in both the momentum balanced 
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(Cho|2005 ) and unbalanced (Cho & Lazarian|2014) situ- 
ations. Even the study of mildly relativistic MHD turbu- 
lence is in its infanc y, having only been treated so far in 
a handful of studies dZhang; et aL|2QQ^ [Inoue et al.|2QlQ 
Zrake & MacFadyen |2Ull[|2UT3{|^rake|2U14[ ). There are 


by comparison, a great number of Mewtoni an MHD tur¬ 
bulence studies (see e.g. Tobias et al.||2011 , for a review) 
treating all different circumstances, including turbulent 
relaxation. Comparisons with them will be made wher¬ 
ever possible. 

One of the issue s we wi l l expl ore in this paper is the 
applicability of the Taylor (1974) hypothesis to magnetic 
relaxation in force-tree electrodynamics. Taylor’s origi¬ 
nal conjecture was that magnetic relaxation would uni¬ 
versally settle in the lowest energy configuration allowed 
by the conservation of total magnetic helicity 


/ 


H = A‘ Bd'^x. 


( 1 ) 


These so-called Taylor states are linear force-free equi¬ 
libria, having electric current density J that is not only 
aligned with the magnetic field, but is also uniformly 
proportional to it, i.e. they solve the constraint 

V X B = aB (2) 

for a global inverse length scale a. Such field config¬ 
urations are monochromatic^ all their magnetic energy 
is concentrated around the spatial frequency a. The 
converse of Taylor’s conjecture is that relaxation may, 
in some circumstances, end in a more general force-free 
equilibrium in which a could vary from one magnetic 
field line to another. In such non-linear equilibria, the 
highest values of a are associated with the smallest scale 
coherent structures, which may be current layers or flux 
tubes, and are associated with peaks in the intensity of 
electrical current flow. 

Counterexamples to Taylor’s conjecture do exist, but 
those identified so far apply to settings in which gas pres¬ 
sure plays a role. For example, hydromag netic relaxed 
states with non- uniform a were repo rted by |Amari fc Lu-| 


ciani 


(2000|) and|Pontin et al.|(|2013|) where the magnetic 
held lines terminate on conducting plates, a boundary 
condition that is motivated by the physics of the solar 
corona. More general hydromagnetic equilibria have also 
been found in simulations of stratified environments such 
as stellar interiors, a setting that has been ext ensively ex¬ 


plored b y iBraithwaitel (|2Q Q6[ |2QQ8[ 
( 2Q1Q| ). iGruzinov "j 2U09[ ) fol lowed 


2008 2009) and Duez et al. 


owed incompressible MHD 
relaxation of a non-helical magnetic field in two dimen¬ 
sions [^and found that it did not decay toward the Taylor 
minimum (total annihilation of the field in this case), but 
instead was halted in an approximate equilibrium with 
many current layers, beyond which further decay was 
only made possible by slow resistive evolution. 

Our study makes frequent use of the periodic short- 
wavelength Taylor states as initial conditions. A Tay¬ 
lor state of frequency and helicity H has an energy 
aoHl2^ a fraction l — ai/a^ of which could be dissipated 

^ Gruzinov’s two-dimensional simulations followed only the in¬ 
plane magnetic field. In the rest of this paper, “two-dimensional” 
means that translational symmetry is enforced along the z-axis, 
but Bz need not vanish. This setting is sometimes referred to as 
2.5D. 


without changing the total helicity (where = 27r/I/ is 
the lowest allowed frequency, although we will use L = 27r 
so that ai = 1). This implies that their free energy 
supply can be arbitrarily large, and so raises the ques¬ 
tion of their mechanical stability. Very recently, [East 
et al. (|2015 ) found that in FFE as well as in relativis- 
tic MHD, generic examples of the 3D, periodic o^o > 1 
Taylor states are unstable to small, ideal perturbations, 
with a growth rate that is proportional to the inverse 
Alfven time. Upon saturation of the linear instability, 
decay enters a turbulent stage that lasts until the re¬ 
maining energy a\H/2 resides at the lowest allowed fre- 
quency a^. This b ehavior bears out the predictions of 


Erisch et al. (|1975) which were based on the prediction 
that turbulence would generically shift magnetic helicity 
toward large scales. 

Conventionally, this so-called inverse cascade has been 
thought to operate efficiently only when the field is 
strongly helical, a belief which has dramatic conse- 


quences fo r large-scale dynamo theory (Blackman & 


fields since the early univer se 


Banerjee & Jedamzik 2004) 


(jOlesen 

1997 

Son 

1999 

iowever, an e; 

tiicient in- 


verse cascade was recently observed to occur even when 
the field was fully no n-helical, in both Newtonian (Bran¬ 
denburg et al.||2015|) and relativistic (|Zrake||2014|) MHD 
settings. Although the magnetic energy eventually de¬ 
cays toward zero, the relaxation evolves in a self-similar 
manner, depositing energy in structures larger than the 
coherence scale which increases over time until it at¬ 
tains the system size. In this study we will show that all 
settings of freely decaying turbulence in force-free elec¬ 
trodynamics, 2D and 3D, helical and non-helical, exhibit 
inverse cascading. The 2D and helical case is particularly 
fast and nearly conservative; as time goes on, magnetic 
energy is shifted toward ever-increasing scales while suf¬ 
fering a diminishing rate of dissipative losses. 

Our paper is organized as follows. We briefly describe 
the theory of force-free electrodynamics and its invari¬ 
ants in Section There we also discuss the special case 
of two dimensions, and define the additional topological 
invariants that it imposes. We then outline the numer¬ 
ical scheme that is used to solve the EEE equations in 
Section and describe our numerical implementation 
of various diagnostics such as power spectra, character¬ 
istic scales, and the helicity invariants. Our simulation 
results, including the energy of relaxed magnetic con¬ 
figurations, an analysis of coherent structures, spectral 
energy distributions, and details of the inverse cascade, 
are presented in Section We discuss the implications 
of these results for astrophysical gamma-ray sources in 
Section!^ and also point out how our results might aid in 
the interpretation of two-dimensional (including axisym- 
metric) calculations. Appendix [A| contains some details 
on the numerical convergence oiour scheme. Through¬ 
out our paper, we use units in which the speed of light 
c = 1. The domain scale L is set to 27r so that the small¬ 
est spatial frequency is 1, and time is reported in units 
of the light-crossing time Ljc. 

2. FORCE-FREE ELECTRODYNAMICS AND ITS 
INVARIANTS 

Eorce-free electrodynamics describes the flow of elec¬ 
tromagnetic energy in a charge-supplied medium with 
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TABLE 1 
Summary of runs 


Grid resolution 


5123 

512^ 

5123 

512^ 

3072^ 

40902,61442^ 81922,122882 
163842 

2562,5122,10242, 20482, 40962 
40902 
122882 
122882 
10243 
10243 

163842 

163842 

7683 

7683 


e 

ao 


See Figure 

1 

16 


[ 


1 

16 



1 

1 


16 


2 

0 


16 


2 

1 

256 

1 



1 

256 



1 l—l 

1 

256 




1 

8,16,32,64,128 



1 

32 




0 

V65641 R 

i 256 


H ^ 

1 

U6564T 

i 256 

2_r 

n[m 

0 

\/2305 

i 48 

n 


1 

\/2305 

i 48 

TT 

fuji 

0 

VesMT R 

i 256 



1 

U6564T 

i 256 


n 

0 

V2^ 

i 48 


n 

1 

U2^ 

i 48 


rr 


Note. — Shown is a summary of the runs along with the figures 
they are referenced by. Helical runs have e = 1 and non-helical 
runs have e = 0. Those whose initial data is the 2D ABC field have 
ao values that are exact integers, whereas randomized initial data 
have ao values that are not exact integers. 


vanishing inertia. This approximation is useful for 
plasma environments in which the energy density of the 
magnetic field greatly exceeds contributions from the 
matter. Under such conditions, the flow of electrical cur¬ 
rent responds rapidly to changes in E and B in order 
to cancel the Lorentz force density pE + J x B, where 
p = V • E is the net electrical charge per unit volume. 
FFE thus admits an Ohm’s law tha t is a function of E 
and B alone, (e.g. McKinn^|2QQ6a] ) 


J = ^(B.VxB-E-VxE) 


E X B 

52 


P, (3) 


which may be coupled to the Maxwell equations 


dtE= VxB-J 
dtB = -V X E 


(4) 


in order to yield a hyperbolic system of partial differential 
equations governing the ev olution of the six component s 
of the electromagnetic field ( Pfeiffer fc MacFadyen|2Q13 ). 

Evolution of ideal MHD systems in general is con- 
strained by three quadratic invariants — the total en¬ 
ergy U, the magnetic helicity i7, and the cr oss helicitv 
W = / V • Bd^x where v is the fluid velocity (Bekenstein 
1987| . As a limiting case of MHD, force-free electrocLy- 


namics shares these invariants, with the exception of the 
cross helicity, since EEE does not define the fluid veloc¬ 
ity in the direction of B. Most relevant to our study 
of magnetic relaxation is conservation of magnetic helic¬ 
ity, which is a topological invariant of the magnetic field 
alone, and thus functions in the same way for MHD as 
it does for EEE. In both MHD and EEE, id is a robust 
invariant, as it is generally found to be conser ved even in 
the p resence of small non-ideal effects (see e.g. Blackman 

Ml. 


2.1. Energy 


Since the Lorentz force density vanishes in EEE, no 
E • J work is done on the charge carriers and the sys¬ 
tem is formally energy conserving. Nevertheless, time- 
dependent solutions may develop regions in which the 
condition E < B is violated, i.e. no frame exists in which 
the electric field vanishes. This indicates a breakdown of 
the ideal force-free assumption, and the Ohm’s law given 
by Equation must be modified. Although non-ideal 
force- free Ohm’s laws have been pr oposed in the litera¬ 
ture ( Gruzinov||2007 Li et al. 2012), it is still common¬ 
place toTvoIvetEemeiLr^yite^^ until such a breakdown 
occurs, and when it does to simply reduce the magnitude 
of E to prevent E > B. The physical motivation for this 
prescription, which we use here, is that energy is being 
radiated away when charges are accelerated to short out 
the electric field and restore E < B. The numerical evo- 
lution scheme is described in more detail in Section [3Tj 
and in Appendix we confirm that it yields numerical 
convergence of the energy dissipation rate. 


2.2. Topology 

Eorce-free electrodynamics shares the same magnetic 
topological invariants as Newtonian and relativistic 
MHD, the lowest order of which is the total magnetic he¬ 
licity given by Equation Its invariance can be seen as 
stemming from the conservation of magnetic flux through 
a closed field loop, which is why it is commonly re¬ 
ferred to as characterizing the interlocking of the mag¬ 
netic field. Although the helicity density A • B depends 
on the choice of gauge, its integral iLy over any volume 
V bounded by magnetic surfaces is well-defined, and also 
invariant under ide al evolutions, such as Earaday’s law , 
when E-B = 0 (e.g. Brandenburg fc Subramanian|2QQ5 ), 


dtHv = -2 / E • Bd^x. (5) 

Jv 

Equation implies that iLy can still be conserved ap¬ 
proximately in the presence of non-ideal processes, as 
long as the volume in which they occur is small. 

In principle, a domain that admits a partitioning by 
magnetic surfaces has an independently conserved helic¬ 
ity associated with each smaller volume. But in prac¬ 
tice, three dimensional field geometries are too complex 
for such a partitioning to be possible. This is what led 
Taylor to conclude that relaxation is only constrained by 
a single topological invariant — the total helicity. The 
story is different in two dimensions since the magnetic 
surfaces are far simpler; their cross-sections are nested 
closed curves in the x — y plane, and the helicity en¬ 
closed by each functions as an independently conserved 
quantity for as long as that surface retains its identity. 
But non-ideal effects, however small, permit the surfaces 
to merge with one another, erasing their identities and 
shuffling up their conserved helicities. Nevertheless, in 
two dimensions we can still construct a robust topologi¬ 
cal invariant that is far stricter than the total helicity. 

We recall that in the presence of ^^-translational sym¬ 
metry, the in-plane magnetic field is tangent to the iso¬ 
contours of the magnetic flux function Az. So each sub¬ 
volume Vi('0) in which Az > ^ is associated with a 
conserved helicity TLi{'ip) = A • Bd^x. To what¬ 

ever extent the helicities of such volumes remain additive 
with respect to reconnections between their bounding 
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Fig. 1. — Top: Two-dimensional turbulent relaxation in force-free electrodynamics at logarithmically spaced times {t = 
0.08,0.32,1.28,5.12). The initial condition is the oq = 256 ABC field with Bi = 1,B2 = l^Bs = 0 and grid resolution 3072^. Shown 
is the out-of-plane magnetic field component scaled linearly between the initial minimum and maximum values. The small red rectangle 
overlying the right-most panel is the region shown amplified in Figure ^ The end-state is not a linear force-free equilibrium. Bot¬ 
tom: Three-dimensional turbulent relaxation under the same conditions except that olq = 16, the grid resolution is 512^, and the times 
t = 0.625,1.0, 3.0,16.0 are chosen to elucidate the sequence of decay epochs. The color mapping accomodates the instantaneous data range, 
as it decreases appreciably throughout the decay. The end-state is a linear force-free equilibrium with a = 1. 


surfaces, the sum must be a robust 

topological invariant. This assumption is justified by 
Equation!^ when the non-ideal regions (where E • B 0) 
are small. Alternatively, consider the helicity (e.g.lBerger 

fcSd|[T^ - 

4^12 = 24>i4>2 (6) 


of two flux tubes that are once-linked and independently 
non-helical. An example is a cylindrical structure whose 
external, toroidal flux 4>2 wraps an interior, poloidal flux 
4>i. Two such structures that were joined to share an 
outer surface would then enclose a poloidal flux 24>i, 
which is wrapped by the same toroidal flux 4>2. So, to 
whatever extent the “joining” is done without destroying 
magnetic flux, the resulting arrangement has a helicity 
27 ^ 12 . 

The presence of additional invariants is expected 

to place a more restrictive lower bound on the magnetic 
energy than would the total helicity H alone. In par¬ 
ticular, given that Taylor states of the same helicity but 
different wavelength will generally span a distinct range 
of 7/) values, there may be no way for one 2D linear equi¬ 
librium to evolve into another of a different wavelength 
while leaving H('0) unchanged. With this in mind, we 
anticipate that fully relaxed 2D states will have an en¬ 
ergy that exceeds the Taylor minimum of aoHl2. We 
will desc ribe the nume rical determination of in Sec- 
tion |3.3[ and in Section [T3| we confirm that it is conserved 
by directly measuring it in our 2D simulations. 

There is also the question of what should happen to 
configurations in which H = 0 but ^ 0, since such 

a configuration could not attain arbitrarily low energy 
while respecting each of the invariants H('0). Although 
behavior like this may be entirely possible, the partic¬ 


ular ini tial conditions used in this study (described in 
Section |3.2| and summarized in Table generally have 
values of |H('0)| that are much smaller than 2Ub/oli^ so 
we have not been able to observe it yet. Indeed, we will 
see in Sect ion [4^ that our 2D states with zero net helicity 
still decay toward very small energy. 


3. METHODS 

We simulate magnetically dominated, relativistic tur¬ 
bulence on a periodic domain of length L = 27r in ei¬ 
ther two or three dimensions, using solutions of the ideal 
force-free electrodynamics equations, given by Equa¬ 
tion and Equation!^ 


3.1. Numerical scheme 

We evolve the EEE system using a fourth-order finite 
difference scheme. We use standard fourth-order differ¬ 
ence operators to evaluate all the field gradients, and 
standard fourth-order Runge-Kutta time-stepping to ad¬ 
vance the solution in time. 

The EEE system requires three vector constraints to 
be maintained: no monopoles V • B = 0, perfect con¬ 
ductivity E • B = 0, and the existence of a frame in 
which E vanishes E < B. The first two constraints 
are formally preserved by EEE, but can be violated nu¬ 
merically at the level of truncation error. Our scheme 
maintains the solenoidal constraint using the hyperbolic 
diverg ence cleaning scheme proposed by [Dedner et al. 


diverg ence cleaning scneme proposed by [Dedner et al 
( 2002^, and late r used in EEE simulations (e.g. |Palen- 
zuela et al.||2010). This amounts to supplementing F'ara- 
day's law with a magnetic monopole current Jb = — VT, 
where the scalar field T evolves according to the damped 
wave equation 5^^ = “V • B — with r being a non¬ 

physical time scale for quenching the magnetic monopole. 
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E B = 0 is maintained exactly by disregarding the part 
of the truncation error that would give rise to a compo¬ 
nent of E in the direction of B. Numerical noise intro¬ 
duced by finite difference operations can lead to unphys¬ 
ical growth of modes whose wavelengths are comparable 
to the numerical grid spacing. Our scheme suppresses 
these unphysical modes using Kreiss-Oliger dissipation, 
a form of low-pass filtering. Each of the procedures just 
mentioned supplements the FEE equations with terms at 
or below the level of the truncation error, so they do not 
modify the formal convergence order of our numerical 
scheme. 

This numerical scheme was used in East et al. (2015), 
and convergence results, as well as comparisons to rel- 
ativistic MHD simulations and analytical methods can 
be found in that reference. It has been im plemented as 
part of the Mara (Zrake & MacFadyen|2011 ) suite of rela¬ 
tivistic turbulence codes, which has many run-time post¬ 
processing capabilities that allow us to perform spectral 
and statistical analysis of the solution at a high cadence 
while minimizing strain on the host architecture’s filesys¬ 
tem. 


3.2. Initial conditions 

We start our simulations with a monochromatic mag¬ 
netic field, where all the power is at a single wavenumber 
magnitude o^o, and with a vanishing electric field. The 
general expression for our initial conditions is 

B(x) = ^ (ao^k + eik x (7) 

|k|=Q:o 

k • 4'k = 0 
4'k = '4'*_k 


where the parameter e is chosen to be either one or zero, 
corresponding to helical or non-helical configurations, re¬ 
spectively. Helical initial configu ratio ns where ao > 1 are 
unstable equilibria (see Section [tT] ), whereas the non- 
helical configurations are out of equilibrium. Some of our 
initial conditions are randomized, having 
where bk is a random unit vector in the plane orthogo¬ 
nal to k and 0k is a random phase. We also make use of 
a spe cial case of Equation known as the “ABC” solu¬ 
tion (Arnold 1965 Dombre et al. 1986). In general this 
is given by 


( Bs cos aoz — B 2 sin (ao^\ 
Bi cos a^x — Bs sin aoz ] , 
B 2 cos aoy — Bi sin a^xJ 


(8) 


which is highly ordered and fully helical, meaning that 
B = 0^0 A (in the Coulomb gauge). In this study we will 
make frequent use of the case with Bi = B 2 = I and 
H 3 = 0 , which we refer to as the 2 D ABC configuration. 

Our results are based on simulations having a range of 
initial frequencies ao and numerical resolutions — which 
we will refer to by the number of grid points in each lin¬ 
ear dimension N. In general, the quality of our results 
improves when we are able to simulate larger values of ao 
with more separation between the initial length scale and 
the domain length scale. However, features (of size 
in our initial condition need to be resolved by a certain 
number of grid points in order to obtain robust solutions. 


In Appendix |a| we show that 32 cells per are suffi¬ 
cient to keep the error in the global helicity conservation 
smaller than 1%. In 2D we will present simulations with 
ao as large as 256, with resolutions up to 16384^. In 3D, 
we will present simulations with ao as large as 48 and 
resolution 1024^. 


3.3. Diagnostics 

We define the power spectral density of the electric, 
magnetic, and helicity fields, respectively, as 

PE{h) = ^ E Eq-EV2, (9) 

* ki<\ci\<ki-\-Aki 

= i E 

* ki<\ci\<ki-\-Aki 

ki<\ci\<ki+Aki 

where Eq, Bq, and Aq are, respectively, the electric field, 
magnetic field, and vector potential Fourier harmonics of 
wavenumber q. We normalize the Fourier harmonics so 
that the volume integrated electric and magnetic field 
energies Ue and Ub^ and the magnetic helicity H are 
given by 


UE = J2PE{ki)Aki, 

i 

UB = Y,PB{h)Aku 

i 

H = J2Pn{ki)Aki. 


We also define the characteristic frequency of each field 
kE, ks^ and kn as 


j _ EiPx{ki)hAki 
^ EiPx{ki)Aki 


( 10 ) 


where X is one of E, H, or H. The most probable 
wavenumber, where Px{k) is maximal, is denoted by kx- 
In two dimensions, we trac k the “helicity mass” func¬ 
tion discussed in Section 

nii^) = J e{A,{x)-^p)A-Bd^x, ( 11 ) 


where 0 is the Heavyside step function. In practice, 
this diagnostic is more easily computed as the “helicity 
density” function which we calculate by binning 

the lattice points according to their value of and 
assigning the weight A • B. We also create the volume 
distribution dV/d^p by binning points according to Az 
with uniform weights, and the helicity distribution over 
volume dn/dV =%/%■ 


4. RESULTS 

Figure shows the evolution of both two and three 
dimensional freely decaying force-free magnetic turbu¬ 
lence. Both of these calculations are initiated in the 
2 D ABC state, but the one on top takes place on a 
two-dimensional domain where translational symmetry 
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Helicity and volume distribution overV^ 



t-tQ 


Fig. 2. — Total magnetic energy as a function of time (since 
the onset of nonlinear evolution to)- The energy of the terminal 
state takes on one of three values. In the case of two-dimensional 
helical evolution, the end-state contains coherent structures and re¬ 
tains 30% of its initial energy, whereas for three-dimensional evo¬ 
lution the end-state is a longest wavelength Taylor state whose 
energy is reduced by a factor of ao from its initial energy. When 
the decay is non-helical magnetic energy decays perpetually to¬ 
ward zero. The randomized initial condition q:o ~ 16 corresponds 
to Qq = 257 and ao ~ 256 corresponds to = 65641. 


is assumed in the 2 : direction, and the bottom one was 
given a low-level white-noise perturbation to break the 
z-symmetry. The left-most image shows the solution 
shortly after saturati on of the linear ins tability that was 
recently observed in East et a]^j2Q15 ), an overview of 
which is provided in Section |4.l[ The difference be¬ 
tween the two runs is visually evident. While the three- 
dimensional solution becomes increasingly smooth at late 
times, the two-dimensional one maintains a network of 
abrupt field reversals. These structures are force-free ro¬ 
tational cmrent layers, and are examine d in depth in 
Section 14.41 As we will see in Section 14.21 the total 


magnetic energy dissipated is dramatically greater in the 
three-dimensional case than in the two-dimensional case. 
Both runs show evidence of the inverse cascade; large- 
scale coherency of the magnetic field must result from dy¬ 
namical transfer of some magnetic energy toward longer 
wavelengths since the initial spectrum is monochromatic 
around k = a^. The inverse cascade will be examined in 
detail in Section lT6l 


4.1. Linear instability of the exeited Taylor states 

Our helical initial conditions are linear force-free equi¬ 
libria. Clearly they are stable when q^o = 1 since such 
states are global energy minima for a given magnetic 
helicity. The question of the ideal stability of the pe¬ 
riodic shorter wavelen gth (aci > 1) Taylor states has a 
conflicted history (e.g. Mc)ffatt||1986 Galloway & Frisch 


1987 


Er-Riani et al. 


solved in East et al. 


__ which has recently been re- 

j2U15 ). In that study, numerical 


simulations of both FT E arid relativistic MHD revealed 
that 0^0 > 1 Taylor states are linearly unstable to ideal 
perturbations. The instability is marked by exponen¬ 
tial growth of the electric field on roughly the Alfven 
wave crossing time of the initial structure size q^q and 
saturates when the medium attains the Alfven speed. 



- 0.06 - 0.04 - 0.02 0.00 0.02 0.04 0.06 


b 

Fig. 3.— Top — The helicity distribution dHl versus the level 
surface value L of the magnetic potential function, shown at 16 
evenly spaced times up to t = 16 in a 2D ABC ru n wi th ao = 32 
and resolution 4096^. As anticipated in Section |2.2[ d'H/d'0 is 
essentially constant in time. Middle — The volume distribution 
dV/dL: indicating the volume between level surfaces at a given 
L- The separatrix surfaces '0 = 0 initially occupy the greatest 
volume, become the current layers, and end up with the smallest 
share of the volume. Bottom — The helicity per volume dl-L/dV. 
Lighter, wider curves indicate earlier time whereas darker, thinner 
curves indicate later times. The distributions are each normalized 
to unity. 

which for the magnetically dominated case is c, implying 
the existence of regions where E ^ B. This instabil¬ 
ity affects generic states, and the only counter-examples 
that were found were one-dimensional ABC states (e.g. 
Bi = 1, ^2 = 0, ^3 =0) that are stable for all values of 
0 ^ 0 . Such states are pure plane waves having circular po¬ 
larization, and are force-free by virtue of having uniform 
magnetic pressure. All of our helical initial conditions are 
short wavelength and either two or three dimensional, 
and turbulent relaxation begins after the saturation of 
the ideal instability. 

4.2. Energy of fully relaxed eonfigurations 

Here we discuss the magnetic energy associated with 
the end-state magnetic configurations. Since the Taylor 
states have B = o^oA, their energy is given simply by 
aoH/2. In other words, their energy is ao times larger 
than the theoretical lower limit imposed by assuming the 
state reaches a = I at constant H. Whether or not dy¬ 
namical relaxation processes settle with the field in this 
global energy minimum remains an open question, but 
here we provide some evidence to support the following 
conjecture: “force-free magnetic relaxation starting from 
periodic Taylor states ends in the lowest energy configu¬ 
ration allowed by helicity conservation if and only if the 
domain is three-dimensional.” 

We have carried out a suite of calculations belonging 
to one of four categories, being either two or three di¬ 
mensional, and either helical or non-helical. Those that 
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Fig. 4.— Amplification of the small red rectangle overlying the rightmost panel of Figure showing relief plots of the x, y, and 2 ; 
components of the magnetic field (on the top three panels). The bottom panel shows one-dimensional profiles, taken along the horizontal 
centerline (dashed magenta). 



Fig. 5.— Probability density function P{o') at logarithmically 
spaced times. The right-most vertical dashed line is the frequency 
ccQ = 256 of the initial field configuration. The local maximum 
at a 128 is the frequency ac of the current sheets (which is 
resolution dependent), and the maximum at a = 0 corresponds to 
zero-current characterizing the flux domain interiors. 

are non-helical are, by construction, out of equilibrium at 
t = 0 and could decay until they reach zero energy since 
helicity conservation does not place any lower bound on 
their magnetic energy. The helical ones are initially at 
an unstable equilibrium, enter a period of turbulent re¬ 
laxation, and settle in a force-free equilibrium of lower 
energy. We considered initial conditions that are both of 
the randomized type given by Equation and the ABC 
type given by Equation 

Our results are summarized in Eigurej^ which shows 
the time series Ub for each of six different runs. As ex¬ 
pected, the non-helical configurations decay toward zero 
energy in both two and three dimensional settings. In 
three dimensions, the helical configurations all terminate 
in the lowest energy state allowed by helicity conserva¬ 
tion]^ Both a randomized setup with = 257 and the 
2D ABC setup with aQ = 256 showed the same general 

^ The sta te of lowest allowed energy can be referred to (e.g|Sliats| 
|et al.|2005| > as the spectral condensate. 


Current layer frequency versus initial frequency and resolution N 




Fig. 6 .— Depenency of the current layer frequency etc on the 
frequency cto of the initial field configuration (top) and the grid 
resolution N (bottom), etc is defined to be the second local maxi¬ 
mum (other than et = 0) in the probability density function P(et), 
where et = B • V x B/P^. 

behavior. The latter setup was intentionally chosen to be 
identical with the two-dimensional setup apart from the 
inclusion of a low-level (one part in 10^) white-noise per¬ 
turbation introduced to break the ^^-translational sym¬ 
metry. 
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Radial profile of magnetic bubble 



Fig. 7.— Example best-fit (shown in red) to magnetic power 
spectrum data (+’s) for a two-dimensional helical model with ao = 
256 and resolution 16384^. Psik) is shown compensated by the 
best-fit spectral index s = 1.94 (top) and uncompensated (bottom). 
Vertical dashed lines indicate the window used for the fit, between 
the peak magnetic frequency ks and the spectral cutoff frequency 
ki. 

Both of the helical two-dimensional runs terminate 
their relaxation with an energy that is decreased to only 
30% of its initial value. For example, a randomized 2D 
initial condition with ~ 256 settles in a state that has 
roughly 77 times more magnetic energy than the Taylor 
minimum energy state. This is not unique, as actually all 
of our helical 2D runs where q^o ^ 1 (including ag = 16) 
settle in a state whose energy is decreased to roughly 
30%. In Appendix we show that, for the 2D ABC 
0^0 = 256 setup, the nnal energy is numerically converg¬ 
ing to a value very near 30%. This means that the termi¬ 
nal states in two dimensions are not linear equilibria. In 
other words, they do approximately solve the force-free 
condition Equation but the value of a may vary from 
one field line to another. 



0 1 2 3 4 5 6 

dimensionless radius f 


Fig. 8 .— Top — The magnetic pressure profile of a magnetic 
bubble in the dimensionless cylindrical radius r = a^r. The blue 
shaded region indicates the azimuthal standard deviation at each 
radius from the center, and the dashed line shows the best-fit model 
parameters f or t he Gaussian mangetic pressure enhancement given 
by Equation |15| Bottom — The azimuthally averaged value of ol 
along with its predicted value (dashed line) given by Equation |16[ 
The top and bottom insets show two-dimensional relief plots of ub 
and o, respectively. 

domain, namely the global maximum of \Az\ which we 
denote by V^max- For two different linear equilibria with 
frequencies ao and a'o to have the same helicity distribu¬ 
tion, it would be necessary for them to at least share the 
same values of H and requiring that 

(^ Bz 

Woj 


4.3. Helicity distribution 

The fact that the 2D configurations do not relax into 
linear equilibria stems from invariance of the hel icity 


distribution which we introduced in Section 2.2 

Figure confirms that it does indeed remain constant 
over time to a very good approximation, even while the 
volume distribution over the magnetic potential dV/d^^ 
changes significantly. The feature evident in the bottom 
panel of Figure where in the relaxed state most of the 
volume is occupied by the extreme values of the magnetic 
potential, can be connected to the formation of current 
layers which we discuss in the next section. 

As an illustration that invariance of the helicity distri¬ 
bution might prevent 2D relaxation from attaining longer 
wavelength linear equilibria, consider that the magnetic 
potential function Az of any 2D linear equilibrium satis¬ 
fies 


A. = 


H 


2aoL^ 


1/2 


Bz 


where Bz is the root mean square value of Bz and H 
is the total helicity. can be characterized by its 


where Bz^max is the global maximum of \Bz\. Note that 
in general, Bz < Bz^ma^- For the particular case of the 
2D ABC state Bi = B 2 = Equation leads to the 
requirement that > (ao/2, and therefore its relaxation 
into a linear state with twice smaller energy is impossi¬ 
ble. In other words, there is no way for the 2D ABC state 
represented in Figurewhose ao =32, to evolve into an¬ 
other linear equilibrium whose frequency is smaller than 
16, while preserving We suspect this argument 

could be generalized further, but for now we leave it as a 
conjecture that the wavelength of a linear 2D equilibrium 
may be uniquely specified by its helicity distribution — 
which if true would render it impossible for one linear 2D 
equilibrium to relax into another of lower energy. 


4.4. The current layers 

During the turbulent relaxation of helical two- 
dimensional configurations, the solution consists of op¬ 
positely signed flux domains separated from one another 
by a network of rotational force-free current layers. The 
flux domains (black and white regions of Figure have 
nearly uniform Bz and are thus relatively current-free. 
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Across the current layers, the magnetic field direction 
rotates through approximately 180°, while its magni¬ 
tude (and thus magnetic pressure) remains fixed. One 
such current layer is shown in Figure where a one¬ 
dimensional profile has been taken along me x-axis, pass¬ 
ing through the layer where it is aligned with the ^-axis. 

It is evident that the current layers have a character¬ 
istic frequency, which we denote by and determine 
empirically as follows. Since the solution is near a force- 
free equilibrium, the current is J ^ aB for some spatially 
dependent frequency 

a = B-VxB/B2. 

We anticipate that the probability density function P{a) 
will have two local maxima — one at a = 0 correspond¬ 
ing to the potential flux domain interiors and the other 
at the frequency o^c, marking the frequency of the cur¬ 
rent layers. Figure confirms this to be the case. It 
shows P{a) at logarithmically spaced times throughout 
the relaxation, and reveals the location of the second 
peak once the solution is sufficiently close to a force-free 
equilibrium. For this particular run, with = 256 and 
N = 3072, the value of ac ~ 128. 

It turns out to be a coincidence that in this case 
ac ~ <ao/2. We have performed a family of calculations, 
varying ao between 8 and 128, and varying N between 
128 and 8192. For each run, we recorded the value of o^c, 
time-averaged over roughly 100 snapshots between time 
t = 10 and t = 16, which was late enough that the second 
peak in P{a) had emerged in each run. There was no sec¬ 
ular evolution of ac. Figure [^reveals that current layers 
become increasingly narrow with higher resolution, but 
also with increasing initial frequency ao. The scaling is 
consistent with the expression 


ae = (13) 

where ki = N/30 is the turbulence cutoff frequency, 
which has been determined W modeling the magnetic 
energy spectrum, (see Figure!^ and is insensitive to ini¬ 
tial conditions, depending only on the numerical scheme 
and grid resolution. We emphasize that Equation [T^ is 


strictly empirical, in that it matches the data shown in 
Figure 

The scaling given by Equation indicates that in the 
infinite resolution limit, the current layers will have zero 
characteristic length. We are not able to say whether the 
scaling could be associated with a physical property of 2D 
EEE at finite Reynolds number, or if it might depend on 
the details of the numerical scheme. This question could 
be resolved by imposing the turbulence cutoff frequency 
ki explicitly, and then varying the numerical resolution. 
In other words, solutions to a resistive EEE system at a 
given conductivity parameter will be necessary. 


4.5. Solitary magnetic bubbles 

The flux domain interiors contain another type of co¬ 
herent structure. They are long-lived bubbles of toroidal 
magnetic field that are confined by ambient magnetic 
pressure. The circular object to the left of the current 
layer in Eigure is an example. They may be referred 
to as flux tubes or magnetic vortices, but we will refer to 
them as magnetic bubbles since we believe t hey a re sim¬ 
ilar objects to those studied by Gruzinov (2010). The 



Fig. 9. — Magnetic energy spectra Pb{P) shown at logarithmi¬ 
cally spaced times for a 2D simulation of freely decaying force- 
free turbulence. Randomized initial conditions with ao ~ 256 and 
N = 16384, where e = 0 (non-helical) on top and e = 1 (helical) 
on bottom. 


bubbles are helical, force-free magnetic field structures, 
having J very well aligned with B. But they are not lin¬ 
ear equilibria; the value of a varies with distance r from 
the axis, and all the bubbles we examined had similar 
radial profiles a{r). The current flow is parallel to the 
background magnetic field near the core, but an equal 
and opposite return current flows through the sheath. So 
they could also be thought of as co-axial electric current 
channels oriented along the symmetry axis. 

Eorce-free configurations with axial and z-symmetry 
satisfy the ordinary differential equations 

B^a PB'^ = 0 (14) 

- B^a + = 0 

where a{r) needs to be specified to make the solution 
unique, as well as the boundary conditions Bz{r) -> B^o 
and B(j){r) ^ 0 as r ^ oo. Requiring the total current 
/ = / B(i)rd(p to be zero only forces B(f){r) to decrease 
faster than I/r. Alternatively, the magnetic pressure 
UB{r) may be specified instead of a{r). One simple as¬ 
sumption, consistent with measurements of the bubbles, 
is that the magnetic pressure enhancement, relative to 
the ambient pressure 5^/2, is Gaussian. This gives a 
pressure profile 


us{f) = tBl(l + fe-^^) (15) 

in the dimensionless radius r = a^r for a bubble of radius 
where / is the magnetic pressure enhancement rel¬ 
ative to background and could be any positive number. 
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Fig. 10. — Total magnetic energy Ub as a function of time (top row) and evolution of the magnetic frequency ks and, for the helical 
runs, the helicity frequency kn (bottom row). One representative run is chosen from each of the four settings — 2D/3D non-helical, and 
2D/3D helical. 


The corresponding solution of Equation 14 


Bz (f) = Boo ^l + /e-^yi-f2) 


which has the profile 


a{r) = ab 


(2 — f^)e 

a/(1 - + /-I 


(16) 


In Figure we show the radial profile of magnetic 
pressure, azimuthally averaged around the bubble’s axis. 
This is the same object depicted to the left of the current 
layer in Figured We also show the expression given in 
Equation \TE\ with its best-fit model parameters. For this 
particular object, the best-fit fractional magnetic energy 
enhancement was / = 0.41 and the best-fit inverse ra¬ 
dius was 0^5 = 162. The lower panel of Figure shows 
the radial profile of a, and also the model predicted value 
given by Equation pT| with the best-fit parameters. 

So far we have notdetected these structures in 3D sim¬ 
ulations, but that does not mean they never happen. It 
is crucial to address whether they are even stable in three 
dimensions, and if they are, then under what conditions 
they may be attractors. 


4.6. The inverse easeade 

The inverse cascade can be characterized by the rate 
with which the magnetic coherence scale (Equa¬ 
tion flO) migrates toward longer wavelengths. We ob¬ 


serve inverse cascading in force-free electrodynamics for 
every setting we have considered, whether the turbulence 
is 2D or 3D, helical or non-helical. Figureshows evolu¬ 
tion of the magnetic energy spectrum P^k) over time, 
for representative helical and non-helical runs in 2D. We 
note how in both cases, the remaining magnetic energy 
resides at an increasing scale over time. We also note 
the spectral energy density at long wavelengths to be an 


increasing function of time, indicating that seleetive de- 
eay of the short wavelength modes alone cannot explain 
growth of the spatial coherency. This further generalizes 
observations mad e of non-helical 3D MHD turbulence, in 
bot h Newtonian (Brandenburg et al.|2015) and relativis¬ 
tic ( Zrake||2Q14 ) settings. 

Figure |l0| shows the evolution of the characteristic 
magnetic frequency ks for one model from each of the 
four categories. For the helical runs we also show the 
evolution of the helicity frequency kn- Except for the 
general feature that over time, both ks and kn move to 
smaller frequency, there is no single decay law that char¬ 
acterizes all these settings. Some runs exhibit power-law 
time dependence of ks, or kn, but not all. Power- 
law dependence is considered evidence of self-similarity 
in the relaxation process, meaning the soluti on evolves 


only by resca ling itself in space and time (e.g. Landau & 


Lifshitz 1987). Self-similar evolution may occur when the 


characteristic scales are all sufficiently smaller than the 
domain size, so that processes around the coherence scale 
are not yet contaminated by the requirement of period¬ 
icity at the domain scale. For this reason, large values 
of 0 ^ 0 , and thus high domain resolution, may be required 
for self-similarity to emerge. For 3D, freely decayin g non¬ 


helical relativistic MHD turbulence, 


Zrake 


(2014) found 


power-law dependence of oc for ao ^ 48 initial 

conditions. For the same conditions in FFE, ks evolves 
with a power-law index of —0.48, as shown in the second 
column of Figure However, the dependence oi ks 
on time is not convincingly power-law, indicating that a 
scale separation of 48 may not be sufficient to yield 
a decisive measurement of the decay index. In the fu¬ 
ture, we plan to follow-up on this with higher resolution 
simulations. 

In freely decaying 2D non-helical force-free turbulence 
with very high resolution (12288^) and ^ 256, clear 
power-law behavior of is observed, with an index of 
—0.38. The helical 2D case is different. As shown in 
Figure 10 (column 3), the energy drops to its terminal 
value of ~ O.SUsit = 0) long before k^^ approaches the 
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itself over all available scales, with the bulk of the energy 
around ks (which decreases over time as discussed in 


frequency k 

Fig. 11.— Power spectra of the electric (to^ and magnetic (bot¬ 
tom) energy densities as given by Equation for one model from 
each of the catagories 2D/3D and helical/non-helical. The solid 
curves are measured from simulations of whose resolution is 16384^ 
in 2D and 1024^ in 3D, while the coinciding dash-dotted curves 
are measured from simulations whose resolution is slightly lower, 
12288^ a nd 7 68^. Dashed lines indicate the best-fit parameters of 
Equation 1 17| An example fit is shown in more detail in Figure 

domain scale. At times later than t = 0.15, the merging 
of flux domains (evident in Figure moves magnetic 
energy to progressively larger scales vmile suffering slower 
and slower dissipative losses. During this epoch, the both 
the helicity and magnetic frequencies decay with a power- 
law index of roughly —0.55. 

In general, freely decaying magnetic turbulence must 
exhibit an inverse cascade when the magnetic helicity 
is near maximal (H ^ UB/kB), if dt H = 0 is to be 


satisfied (Frisch et al. 1975 Christensson et al. 2001 


Cho| |2011 r However, the same conclusion cannot be 
drawn from magnetic helicity conservation alone when 
H ^ Us/kB^ as is the case for our non-helical runs. 
Observation of the non-helical inverse cascade was thus 


Section 4.6) and a power-law tail extending to the spec¬ 
tral cutoft frequency ki which lies consistently near A^/30 
when the grid resolution is N. We have determined the 
index s of the power-law tail by fitting the iogarithm of 
the spectrai distributions PE{k) and Psik) to the modei 
function 

f{x) = Ae-<-^-^°^ +XS (17) 

where x = log/c, and A, xq, cr, and s are modei param¬ 
eters. We have found it expedient not to try and modei 
the high-frequency cutoffs, we just use frequency bins be¬ 
tween the peak frequency ks and the cutoff ki to obtain 
our fit. A representative spectrai fit is shown in Figure 
Figure [IT] shows PE{k) and Psik) for one modei from 
each of tHe categories 2D/3D and heiicai/non-heiicai 
aiong with the best-fit modei given by Equation We 
use randomized initiai conditions for each modei, and 
the resoiution is 16384^ in 2D and 1024^ in 3D (coincid¬ 
ing dash-dotted curves are taken from simuiations whose 
resoiution is 12288^ and 768^ to indicate numericai con¬ 
sistency). The spectra represent the soiution at an inter¬ 
mediate time, when ~ 8, and we find the power-law 
indices to be quite stable during a window of time be¬ 
ginning on the snapshot chosen for each model. The 
power-law index of the electric field spectral energy dis¬ 
tribution Pe(/c) is between 0.96 and 1.18 for the various 
models, with both of the helical models having an in¬ 
dex of 1.18. The power-law index of the magnetic field 
spectral energy distribution Psik) is different in two and 
three dimensions. In both helical and non-helical 3D set¬ 
tings, it is notably close to the Kolmogorov value of 5/3 
Note that ~ ~ 


Zrake 


(2014) and Brandenburg et al. (2015) 
both measured an index near 2 tor freely decaying non¬ 
helical MHD turbulence. The helical 2D model has an 
index of 1.93, and it is tempting to conclude that the 
true value is 2, but actually the best-fit parameter, for 
aii resoiutions up to 16384^ is significantiy smaiier than 
2, aiways being between 1.92 and 1.96. The non-heiicai 
2D modei is found to have an index of 1.41. 

5. DISCUSSION AND CONCLUSIONS 

Using high resoiution simuiations of the force-free eiec- 
trodynamics equations, we have studied freeiy decaying, 
magneticaiiy dominated reiativistic turbuience for the 
first time. We focused on various differences between 
the magnetic reiaxation process in settings where the do¬ 
main is two and three dimensionai, and where the mag¬ 
netic fieid is heiicai and non-heiicai. We found that he- 
licai, two-dimensionai reiaxation terminates in a state 


seen as a surprise (Zrake|2014 

Brandenburg et al.|2015 

). posed by helicity conservation, and whose volume is pre- 

It was suggested in Zrake (|2U; 
verse cascade may reflect the t< 

L4p that the non-helical in- dominantly current-free, but is punctuated by coherent 

endency for aligned current structures — namely current layers and solitary magnetic 


struct ures to attract one another. Brandenburg et ai 


(2015) found that net transfer of energy from smaii to 


iarge scaies was about twice iarger in heiicai than non- 
heiicai freeiy decaying Newtonian MHD turbuience. 


bubbies. We tried to determine what sets the width of 
the current iayers, and determined that it depends not 
only on the turbulence cutoff frequency (or grid resolu¬ 
tion), but also on the frequency of the initial condi- 


4.7. Power spectrum of electric and magnetic energy 

In all of our initial conditions the magnetic energy 
is concentrated around a single frequency, Psik) oc 
6{k — ao). As turbulence develops, energy redistributes 


^ Between t = 0.6 and t = 1.2, the spectral indices show negli¬ 
gible secular evolution, and a standard deviation with respect to 
different time levels that is below the level of 1%. The spectral 
indices reported are the instantaneous values, which are within a 
standard deviation of the mean. 
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tion. The solitary magnetic bubbles are axisymmetric, 
non-linear force-free equilibria, and are consistent with 
Gaussian magnetic pressure enhancement relative to the 
surrounding relatively current-free volume. The three- 
dimensional stability of these structures remains an open 
question. 

The unusual behavior of 2D relaxations can be under¬ 
stood in terms of additional topological constraints that 
are imposed by the extra symmetry. We proposed that 
unlike generic 3D relaxation, whose topology is only con¬ 
strained by the total helicity invariant i7, 2D relaxations 
are subject to a whole spectrum of helicity invariants 
H('0), one associated with each value of the magnetic 
potential t/;. Although this invariance is only guaranteed 
when magnetic reconnections are restricted to regions 
of zero volume, our simulation data showed that 
remains unchanged throughout the evolution to a very 
good approximation. 

All of the settings we considered exhibited inverse cas¬ 
cading, in which some magnetic energy is redistributed 
toward progressively longer wavelengths. The rate of in¬ 
verse cascading was characterized by the time evolution 
of the magnetic frequency which was found to de¬ 
crease faster when the field is helical than non-helical, 
and also faster in 3D than in 2D. The inverse cascade 
of 2D helical turbulence is nearly conservative; merging 
of magnetic flux domains moves energy to larger scales 
while suffering a diminishing rate of dissipative losses. 
The non-helical inverse cascade has ks ^ ^-o.ss 2D 
an d ks c x in 3D, in rough agreement with results 

of Zrake (2014|) for 3D turbulent relaxation in relativis¬ 
tic aTHD. Better scale separation (larger ao) and thus 
higher numerical resolution is needed to confirm the three 
dimensional scaling measurement. 

5.1. Astrophysical gamma-ray sources and 
magnetoluminescence 

We have examined turbulent relaxation in force-free 
electrodynamics with the motivation of elucidating the 
physical mechanism behind extremely fast time vari¬ 
ability that is characteristic of astrophysical gamma-ray 
emitters, including the Crab Nebula, many blazars, and 
nearly all gamma-ray bursts. The extreme energetics and 
temporal intermittency of gamma radiation from these 
sources require a mechanism in which plasma promptly 
converts the majority of its magnetic energy into high 
energy particles and radiation. Furthermore, the emit¬ 
ting regions are thought to be strongly magnetized, and 
known to lie a great distance from the primary mover 
(pulsar, progenitor star, or black hole). These facts are 
highly suggestive that electromagnetic outflows may con- 


prompt relaxation into the Taylor minimum energy state, 
supporting the idea that magnetic energy can be dissi¬ 
pated completely in a light-travel time. But the same re¬ 
sult suggests that persistent, meta-stable structures are 
not a generic outcome of turbulence in force-free elec¬ 
trodynamics. This is not to say that such behavior is 
impossible, as we have onl y considered a smal l class of 
initial conditions. In fact, Smiet et al. (2015) very re¬ 
cently identified magnetic arrangements that on 3D pe¬ 
riodic domains, in full MHD, relax to non-linear hydro- 
magnetic equilibria. So (1) is possible, at least in the 
hydromagnetic case. The volatility of such objects, and 
the generality of conditions under which they may arise 
remain important questions for the future. 

5.2. Comparison with other studies of magnetic 
relaxation 

In this study we have begun to address the question 
of whether relativistic, force-free magnetic relaxation on 
periodic domains generically ends in a Taylor state, and 
found evidence to support the view that it does, provided 
the domain is three dimensional. But thus far we have 
only considered a restricted class of initial conditions — 
namely isotropic, monochromatic fields that are either 
linear force-free equilibria (helical) or completely non¬ 
helical. 

Several studies in full MHD have now identified set¬ 
tings that relax to more general hydromagnetic equi¬ 
libria for which J x B = Vp, where the Lorentz force 
density balances the gradient of gas pressure p. Ex- 
amples include stratified three-dimensional environments 
(Braithwaite 2006 2008) where the field is helical, and 
two-dimensional periodic sett ings where the fi eld is in¬ 
compres s ible a nd no n-helical ( Gruzinov ||2009 ). Amari & 
Luciani (2000) and Braithwaite (|20I5|) have both pro¬ 
vided examples of 3D hydromagnetic relaxation which 
first develop current layers, and then proceed to a smooth 
configuration via resistive processes. Both studies used 
boundary conditions where at least one of the dir e ctions 
was not periodic. Very recently, Smiet et al. (2015) 
has found instances of 3D hydromagnetic relaxation end- 
ing in smooth, non-linear equilibria with non-uniform 
pressure, even when the boundary is periodic or open. 
Such boundary conditions are highly relevan t fo r the 
astrophysical processes mentioned in Section |5.1[ and 
an important question is whether their results possess 
any force-free analogues. In particular, if stable and 
non-linear force-free equilibria do exist in 3D away from 
boundaries, then how likely are they to arise under real¬ 
istic astrophysical conditions? 


tain persistent magnetic structures with copious free en¬ 
ergy supplies, whose spontaneous disruption could be 
linked to the observed flaring events. A scenario like 
this was referred to in Blandford et al.| (|2014|) as mag- 
netoluminescence. For it to be plausible, it is necessary 
that (1) meta-stable, force-free (or hydromagnetic) equi¬ 

5.3. Effects of imposing extra symmetries in 
astrophysical simulations 

Many studies of relativistic plasma and MHD pro¬ 
cesses are carried out assuming either translational or 
rotational symmetry, including simulations of force-free 

libria can exist far from any supporting boundaries, (2) 

electrodynamics ( 

McKinney ||2006bl iTchekhovskoy 

et al. 

that such objects can form under realistic astrophysical 

2008), relativistic IVIHD (|Barkov &; Komissarov 

:^008 

conditions, and (3) that upon their disruption, magnetic 

Komissarov et al.||2009l |Komissarov V Barkov 

TO 

energy is promptly and completely dissipated. 

Mizuno et al.| 20111), and particle-in-cell (PIC) Simula- 

In this paper, we have begun to address the points 

tions ( bpitkovsky||2008 Keshet et al.||2009). As we have 

(1) and (3) and found results that are at least partially 

seen, 2D magnetic relaxations are far more topologically 


encouraging. All of our periodic 3D simulations exhibit constrained, and persist in configurations of much higher 
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energy than equivalent 3D relaxations. Furthermore, ax- 
isymmetric calculations are expected to exhibit similar 
artificialities to the slab-symmetric ones studied here, 
since they both share the same topological simplifica¬ 
tions. In particular, the axisymmetric magnetic surfaces, 
now toroidal shells that are labeled by their value of the 
azimuthal vector potential component ^4^, each enclose 
a conserved magnetic helicity. 

Our results indicate that as a field’s complexity in¬ 
creases, so does the discrepancy between the energy of its 
most relaxed state in 2D and 3D. So, axisymmetry may 
be appropriate when the field is near a stable force-free 
equilibrium, when little or no energy resides in higher- 
order radial or angular modes. But when these modes 
are populated, the imposed symmetry could make them 
artificially persistent. 

Numerical studies of pulsar magnetospheres and mag- 
netar flares are quite challenging, and much of the 
progress in this field has been obtained by restricting to 
axisymmetry. S i mulations using both FFE dParfr ey et al. 
2Q12| r^|2Q 12| |Parfrey et al. [2013|) and Pi(J (jCerutti 
et al.| |2U15|) show the development of complex struc¬ 
ture m the meridional plane. Similarly, short wave¬ 
length toroidal magnetic structure has been found in 
axisymmetric FFE simula tions of black hole a ccretion. 


Eor example the results of Parfrey et al. (2014) suggest 
that even disordered (as opposed to large-scale) magnetic 
fields advected inward by an accretion disk could facil¬ 
itate angular momentum extraction from a black hole. 
The differences in 2D versus 3D found here suggest that 
extending these studies to be fully three dimensional, as 
that becomes computationally feasible, will be an excit¬ 
ing frontier, and will likely reveal qualitatively new dy¬ 
namics. 

Another setting where 2D and 3D calculations are ex¬ 
pected to differ is shock-generated turbulence, which may 
account for the relatively high magnetizations inferred 
from non-thermal emission spectra of astrophysical shock 
fronts — both in non-relativistic settings such as super¬ 
nova remnants and relativistic settings such as gamma- 
ray burst afterglows. Amplification of the magnetic field 
by turbulence in and around the shock, and its subse¬ 
quent decay in the po st-shock flow has b een studied ex- 
tensively in bot h two (Mizuno et al. |2Q11[ ) and three dlir 


one et al. 2010) dimensions with re lativistic MHD, and 


also in 2D and 3D PIC simulations (Sironi & Spitkovsky 
2009 Sironi et"ar]|2013 ). In this study we have observed 


power-law decay of the magnetic energy in all settings, 
but with a steeper index in 3D than in 2D. This should 
be kept in mind, as even minor differences in the decay 


Eermi processes (iLemoine et al.| 

20061 INiemiec & Os- 

trowski 2006 INiemiec et al. 2006 

■Pel/etier et al. |:iOuD| 

as well as interpretations of GRB 

afterglows ( 

Gruzinov 

& Waxman||1999 Rossi & Rees||2003 Lemoine 

Mir 


also on Comet at the San Diego Supercomputer Cen¬ 
ter (SDSC) through XSEDE grant AST150038, as well 
as Pleiades of the NASA High-End Computing (HEC) 
Program through the NASA Advanced Supercomputing 
(NAS) Division at Ames Research Center. 
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Simulations were run on the Bullet Cluster at SLAC 
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APPENDIX 

NUMERICAL CONVERGENCE 

Here we demonstrate some numerical convergence 
properties of our scheme. We have chosen the con¬ 
servation of magnetic helicity AH{t) and the time se¬ 
ries of magnetic energy as diagnostics. Conver¬ 

gence properties are reported for the 2D helical runs 
with 0^0 = 256 and grid resolutions 4096^, 6144^, 8192^, 
12288^ and 16384^. This configuration was found to be 
representative of convergence properties in other settings 
reported in this paper. All runs were evolved for at least 
two light-crossing times. In order to establish the numer¬ 
ical convergence order n, we model the error of the nu¬ 
merical solution Hh with grid spacing h as = VoEk^ 
where ^ is a constant and yo is the extrapolated solu¬ 
tion. This is similar to a Richardson extrapolation, but 


instead of fitting for the coefficient for each integer 
power of h, we have fit for the single error coefficient E 
and convergence order n. 

The upper left panel of Figure [T^ shows the fractional 
change in magnetic helicity H{t)jHo — 1 as a function 
of time for each resolution. For resolutions > 4096^, the 
helicity change is never worse than ±10%. For > 8192^ 
it is never worse than ±1%, and for > 12288^ it is never 
worse than ±0.1%. The extrapolated value of the helic- 
ity change is consistently a gain of about 0.1%, and the 
con verg ence order (shown on the lower left panel of Fig¬ 


ure 12) is between 2.8 and 2.9. The right panel shows 
the evolution of magnetic energy Usit) at each resolu¬ 
tion. Less dissipation occurs for each higher resolution, 
but the sequence converges consistently at first order. 
The extrapolated value of the magnetic energy at t = 2 
is 0.30, and it changes by less than 1% between t = 1 
and t = 2. 
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Numerical convergence of helicity conservation Numerical convergence of total magnetic energy 



Fig. 12.— Left - Error in the conservation of magnetic helicity H. The upper panel shows the fractional helicity change Ah(t) = 
— 1 on symmetric logarithmic axes (to account for anomolous helicity change of either sign) for six different values of the mesh 
spacing h. The dashed magneta line shows the Richardson-extrapolated value of which remains constant at roughly 10“^. The 

lower panel shows the convergence order of Ah{t) at representative times. Right - Evolution of the total magnetic energy Ub (t) for the 
same six values of the mesh spacing. The dashed magenta line on the upper panel shows the Richardson-extrapolated time series of Usit)-, 
and the convergence order is shown on the lower panel. 
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